Lake Elsinore image processing

library(knitr)
opts_chunk$set(warning = FALSE, message = FALSE, fig.path = 'figs/', dev.args = list(bg = 'transparent'), eval = T)

library(raster)
library(sf)
library(tidyverse)
library(mapview)
library(proj4shortcut)
library(plotly)

prj_geo <- geo_wgs84
prj_pro <- utm_wgs84('11s')

August 21st

Map of in situ samples with chlorophyll concentration.

# get coverage of in situ points by date
insit <- read_csv('raw/CHL_2017_08_21.csv', col_names = F) %>% 
  rename(
    site = X1, 
    lat = X2, 
    lon = X3, 
    chl = X4
  ) %>% 
  mutate(lat = ifelse(site == 34, 33.66954, lat)) 
coordinates(insit) <- ~ lon + lat
proj4string(insit) <- prj_geo
mapview(insit, zcol = 'chl', legend = T) 

Overlap of in situ samples (red numbers) with the spatial extent of each image. Note that the spatial extent is rectangular and not all pixels in the extent have values.

# transform insit to projected
insit <- spTransform(insit, CRS(prj_pro)) 

# raster tiffs
imgs_pth <- list.files('img/GeoTiff_08_21/', pattern = '*\\.tif$', full.names = T)

# get raster extents
out <- list()
cmbs <- for(i in 1:length(imgs_pth)){
  
  # raster aggregate, dissolve to polygon
  rst_chk <- stack(imgs_pth[[i]]) %>% 
    .[[1]] %>% 
    aggregate(fact = 10, fun = mean)
  values(rst_chk) <- ifelse(is.na(values(rst_chk)), NA, 1)
  ply <- rst_chk %>% 
    rasterToPolygons(dissolve = T)
  
  out[[i]] <- ply
  
}

# fortify polygons for plot
allply <- out %>% 
  map(., function(x){
    dat <- fortify(x) %>% 
      mutate(fl = gsub('\\.[0-9]$', '', names(x@data)))
    return(dat)
    }) %>% 
  enframe %>% 
  unnest

# plot
p <- ggplot() + 
  geom_polygon(data = allply, aes(x = long, y = lat, group = fl), alpha = 0.5) + 
  geom_text(data = data.frame(insit), aes(x = lon, y = lat, label = site)) + 
  coord_equal()
ggplotly(p)

Extracted pixel values for each of four bands in each image where sample sites occurred within an image.

# extract values from raster
out <- list()
for(i in 1:length(imgs_pth)){
  
  # cat(i, 'of', length(imgs_pth), '\n')
  
  # extract raster cells by points in insit
  rst <- stack(imgs_pth[[i]])
  crs(rst) <- prj_pro
  rst_chk  <- rst %>% 
    raster::extract(insit, buffer = 0.5) %>% 
    enframe %>% 
    bind_cols(insit@data, .) %>% 
    select(-name) %>% 
    filter(map(.$value, is.matrix) %>% unlist) %>% 
    mutate(value = map(value, function(x){

      sumpt <- x %>% 
        as.data.frame %>% 
        gather('img_bnd', 'din', everything()) %>% 
        group_by(img_bnd) %>% 
        summarise(
          ncell = length(img_bnd),
          din = mean(din)
        )
      return(sumpt)
      
    })) %>% 
    unnest

  # append to output
  out[[i]] <- rst_chk

}

exts <- out %>% 
  do.call('rbind', .) %>% 
  remove_rownames() %>% 
  separate(img_bnd, c('img', 'bnd'), sep = '\\.') %>% 
  filter(!is.na(din)) # check this
exts

## # A tibble: 240 x 6
##     site   chl img                    bnd   ncell   din
##    <int> <dbl> <chr>                  <chr> <int> <dbl>
##  1    31   189 IMG_170821_213834_0008 1        13 0.314
##  2    31   189 IMG_170821_213834_0008 2        13 0.166
##  3    31   189 IMG_170821_213834_0008 3        13 0.142
##  4    31   189 IMG_170821_213834_0008 4        13 0.141
##  5    31   189 IMG_170821_213840_0009 1        16 0.328
##  6    31   189 IMG_170821_213840_0009 2        16 0.175
##  7    31   189 IMG_170821_213840_0009 3        16 0.145
##  8    31   189 IMG_170821_213840_0009 4        16 0.144
##  9    31   189 IMG_170821_213846_0010 1        17 0.328
## 10    31   189 IMG_170821_213846_0010 2        17 0.170
## # ... with 230 more rows

Scatterplots of band values (din) with measured chlorophyll.

toplo <- exts %>% 
  group_by(site, chl, bnd) %>% 
  summarise(din = mean(din))
ggplot(toplo, aes(x = din, y = chl)) + 
  geom_point() + 
  facet_wrap(~bnd, ncol = 2) + 
  stat_smooth(method = 'lm') + 
  theme_bw()

Working with model selection:

library(caret)
library(MuMIn)

tomod <- exts %>% 
  mutate(bnd = paste0('b', bnd)) %>% 
  spread(bnd, din)

glb <- lm(chl ~ (b1 + b2 + b2 + b3 + b4)^3 + I(b1^2) + I(b2^2) + I(b3^2) + I(b4^2), tomod, 
          na.action = 'na.pass')

tmp <- dredge(glb, evaluate = F)

# createFolds(1:nrow(tomod), k = 5)

September 6th

Map of in situ samples with chlorophyll concentration.

# get coverage of in situ points by date
insit <- read_csv('raw/CHL_2017_09_06.csv', col_names = F) %>% 
  rename(
    site = X1, 
    lat = X2, 
    lon = X3, 
    chl = X4
  ) %>% 
  mutate(lon = ifelse(site == 39, -117.3597, lon))
coordinates(insit) <- ~ lon + lat
proj4string(insit) <- prj_geo
mapview(insit, zcol = 'chl', legend = T) 

Overlap of in situ samples (red numbers) with the spatial extent of each image. Note that the spatial extent is rectangular and not all pixels in the extent have values.

# transform insit to projected
insit <- spTransform(insit, CRS(prj_pro)) 

# raster tiffs
imgs_pth <- list.files('img/GeoTiff_09_06/', pattern = '*\\.tif$', full.names = T)

# get raster extents
out <- list()
cmbs <- for(i in 1:length(imgs_pth)){
  
  # raster aggregate, dissolve to polygon
  rst_chk <- stack(imgs_pth[[i]]) %>% 
    .[[1]] %>% 
    aggregate(fact = 10, fun = mean)
  values(rst_chk) <- ifelse(is.na(values(rst_chk)), NA, 1)
  ply <- rst_chk %>% 
    rasterToPolygons(dissolve = T)
  
  out[[i]] <- ply
  
}

# fortify polygons for plot
allply <- out %>% 
  map(., function(x){
    dat <- fortify(x) %>% 
      mutate(fl = gsub('\\.[0-9]$', '', names(x@data)))
    return(dat)
    }) %>% 
  enframe %>% 
  unnest

# plot
p <- ggplot() + 
  geom_polygon(data = allply, aes(x = long, y = lat, group = fl), alpha = 0.5) + 
  geom_text(data = data.frame(insit), aes(x = lon, y = lat, label = site)) + 
  coord_equal()
ggplotly(p)

Extracted pixel values for each of four bands in each image where sample sites occurred within an image.

# extract values from raster
out <- list()
for(i in 1:length(imgs_pth)){
  
  # cat(i, 'of', length(imgs_pth), '\n')
  
  rst <- stack(imgs_pth[[i]])
  crs(rst) <- prj_pro
  rst_chk  <- rst %>% 
    raster::extract(insit) %>% 
    data.frame(insit@data, .) %>%
    gather('img_bnd', 'din', -site, -chl) %>% 
    na.omit

  # out <- c(out, list(tmp))
  out[[i]] <- rst_chk

}

exts <- out %>% 
  do.call('rbind', .) %>% 
  remove_rownames() %>% 
  separate(img_bnd, c('img', 'bnd'), sep = '\\.')
exts

##     site    chl                    img bnd       din
## 1      1  62.62 IMG_170821_213858_0012   1 0.3521872
## 2     16 380.14 IMG_170821_213858_0012   1 0.3669254
## 3     17 295.34 IMG_170821_213858_0012   1 0.2611819
## 4      1  62.62 IMG_170821_213858_0012   2 0.1900249
## 5     16 380.14 IMG_170821_213858_0012   2 0.1846308
## 6     17 295.34 IMG_170821_213858_0012   2 0.1613828
## 7      1  62.62 IMG_170821_213858_0012   3 0.1485732
## 8     16 380.14 IMG_170821_213858_0012   3 0.1278632
## 9     17 295.34 IMG_170821_213858_0012   3 0.1200351
## 10     1  62.62 IMG_170821_213858_0012   4 0.1408067
## 11    16 380.14 IMG_170821_213858_0012   4 0.1263729
## 12    17 295.34 IMG_170821_213858_0012   4 0.1322991
## 13     1  62.62 IMG_170906_182425_0023   1 0.3257019
## 14    16 380.14 IMG_170906_182425_0023   1 0.3136333
## 15     1  62.62 IMG_170906_182425_0023   2 0.1682546
## 16    16 380.14 IMG_170906_182425_0023   2 0.1731694
## 17     1  62.62 IMG_170906_182425_0023   3 0.2130444
## 18    16 380.14 IMG_170906_182425_0023   3 0.2025952
## 19     1  62.62 IMG_170906_182425_0023   4 0.1939986
## 20    16 380.14 IMG_170906_182425_0023   4 0.1948591
## 21     1  62.62 IMG_170906_182431_0024   1 0.3138918
## 22    16 380.14 IMG_170906_182431_0024   1 0.3464889
## 23     1  62.62 IMG_170906_182431_0024   2 0.1581172
## 24    16 380.14 IMG_170906_182431_0024   2 0.1802249
## 25     1  62.62 IMG_170906_182431_0024   3 0.2159136
## 26    16 380.14 IMG_170906_182431_0024   3 0.2001033
## 27     1  62.62 IMG_170906_182431_0024   4 0.2123534
## 28    16 380.14 IMG_170906_182431_0024   4 0.1866015
## 29     1  62.62 IMG_170906_182437_0025   1 0.3141099
## 30    16 380.14 IMG_170906_182437_0025   1 0.2983865
## 31    31 254.50 IMG_170906_182437_0025   1 0.3381861
## 32     1  62.62 IMG_170906_182437_0025   2 0.1613249
## 33    16 380.14 IMG_170906_182437_0025   2 0.1533507
## 34    31 254.50 IMG_170906_182437_0025   2 0.1642923
## 35     1  62.62 IMG_170906_182437_0025   3 0.2072103
## 36    16 380.14 IMG_170906_182437_0025   3 0.1837897
## 37    31 254.50 IMG_170906_182437_0025   3 0.2112308
## 38     1  62.62 IMG_170906_182437_0025   4 0.1993446
## 39    16 380.14 IMG_170906_182437_0025   4 0.1752855
## 40    31 254.50 IMG_170906_182437_0025   4 0.2066646
## 41     1  62.62 IMG_170906_182443_0026   1 0.3077186
## 42    16 380.14 IMG_170906_182443_0026   1 0.3145586
## 43    31 254.50 IMG_170906_182443_0026   1 0.3824544
## 44     1  62.62 IMG_170906_182443_0026   2 0.1680164
## 45    16 380.14 IMG_170906_182443_0026   2 0.1602393
## 46    31 254.50 IMG_170906_182443_0026   2 0.1722880
## 47     1  62.62 IMG_170906_182443_0026   3 0.2157355
## 48    16 380.14 IMG_170906_182443_0026   3 0.2072434
## 49    31 254.50 IMG_170906_182443_0026   3 0.2396860
## 50     1  62.62 IMG_170906_182443_0026   4 0.2234698
## 51    16 380.14 IMG_170906_182443_0026   4 0.1862677
## 52    31 254.50 IMG_170906_182443_0026   4 0.2482327
## 53     1  62.62 IMG_170906_182449_0027   1 0.3003604
## 54    16 380.14 IMG_170906_182449_0027   1 0.3132594
## 55    17 295.34 IMG_170906_182449_0027   1 0.3272939
## 56    31 254.50 IMG_170906_182449_0027   1 0.4016975
## 57     1  62.62 IMG_170906_182449_0027   2 0.1576792
## 58    16 380.14 IMG_170906_182449_0027   2 0.1483788
## 59    17 295.34 IMG_170906_182449_0027   2 0.1645345
## 60    31 254.50 IMG_170906_182449_0027   2 0.1726661
## 61     1  62.62 IMG_170906_182449_0027   3 0.2268322
## 62    16 380.14 IMG_170906_182449_0027   3 0.2526465
## 63    17 295.34 IMG_170906_182449_0027   3 0.2137336
## 64    31 254.50 IMG_170906_182449_0027   3 0.2233939
## 65     1  62.62 IMG_170906_182449_0027   4 0.1983301
## 66    16 380.14 IMG_170906_182449_0027   4 0.2199366
## 67    17 295.34 IMG_170906_182449_0027   4 0.2050995
## 68    31 254.50 IMG_170906_182449_0027   4 0.2229731
## 69    17 295.34 IMG_170906_182455_0028   1 0.2903806
## 70    31 254.50 IMG_170906_182455_0028   1 0.4145821
## 71    32 214.83 IMG_170906_182455_0028   1 0.0000000
## 72    17 295.34 IMG_170906_182455_0028   2 0.1540359
## 73    31 254.50 IMG_170906_182455_0028   2 0.1831722
## 74    32 214.83 IMG_170906_182455_0028   2 0.0000000
## 75    17 295.34 IMG_170906_182455_0028   3 0.1966201
## 76    31 254.50 IMG_170906_182455_0028   3 0.2158420
## 77    32 214.83 IMG_170906_182455_0028   3 0.2106364
## 78    17 295.34 IMG_170906_182455_0028   4 0.2066191
## 79    31 254.50 IMG_170906_182455_0028   4 0.2331240
## 80    32 214.83 IMG_170906_182455_0028   4 0.1957703
## 81     2 486.96 IMG_170906_182501_0029   1 0.3103488
## 82    17 295.34 IMG_170906_182501_0029   1 0.2969343
## 83    32 214.83 IMG_170906_182501_0029   1 0.4165199
## 84     2 486.96 IMG_170906_182501_0029   2 0.1686566
## 85    17 295.34 IMG_170906_182501_0029   2 0.1491194
## 86    32 214.83 IMG_170906_182501_0029   2 0.1849540
## 87     2 486.96 IMG_170906_182501_0029   3 0.2082499
## 88    17 295.34 IMG_170906_182501_0029   3 0.2056280
## 89    32 214.83 IMG_170906_182501_0029   3 0.2059632
## 90     2 486.96 IMG_170906_182501_0029   4 0.2313217
## 91    17 295.34 IMG_170906_182501_0029   4 0.2107968
## 92    32 214.83 IMG_170906_182501_0029   4 0.2166646
## 93     2 486.96 IMG_170906_182507_0030   1 0.3085389
## 94    17 295.34 IMG_170906_182507_0030   1 0.3139217
## 95    18 305.90 IMG_170906_182507_0030   1 0.2996119
## 96    32 214.83 IMG_170906_182507_0030   1 0.4486934
## 97     2 486.96 IMG_170906_182507_0030   2 0.1510785
## 98    17 295.34 IMG_170906_182507_0030   2 0.1526142
## 99    18 305.90 IMG_170906_182507_0030   2 0.1638244
## 100   32 214.83 IMG_170906_182507_0030   2 0.1995005
## 101    2 486.96 IMG_170906_182507_0030   3 0.1965176
## 102   17 295.34 IMG_170906_182507_0030   3 0.2199249
## 103   18 305.90 IMG_170906_182507_0030   3 0.2149738
## 104   32 214.83 IMG_170906_182507_0030   3 0.2212107
## 105    2 486.96 IMG_170906_182507_0030   4 0.2272809
## 106   17 295.34 IMG_170906_182507_0030   4 0.2050937
## 107   18 305.90 IMG_170906_182507_0030   4 0.1843024
## 108   32 214.83 IMG_170906_182507_0030   4 0.2479621
## 109    2 486.96 IMG_170906_182513_0031   1 0.3167584
## 110   18 305.90 IMG_170906_182513_0031   1 0.3033060
## 111   32 214.83 IMG_170906_182513_0031   1 0.4871817
## 112    2 486.96 IMG_170906_182513_0031   2 0.1560977
## 113   18 305.90 IMG_170906_182513_0031   2 0.1605711
## 114   32 214.83 IMG_170906_182513_0031   2 0.1999720
## 115    2 486.96 IMG_170906_182513_0031   3 0.2126312
## 116   18 305.90 IMG_170906_182513_0031   3 0.2118794
## 117   32 214.83 IMG_170906_182513_0031   3 0.2445490
## 118    2 486.96 IMG_170906_182513_0031   4 0.2171975
## 119   18 305.90 IMG_170906_182513_0031   4 0.2094065
## 120   32 214.83 IMG_170906_182513_0031   4 0.2229194
## 121    2 486.96 IMG_170906_182519_0032   1 0.3505707
## 122   18 305.90 IMG_170906_182519_0032   1 0.3049344
## 123   33 151.64 IMG_170906_182519_0032   1 0.4672850
## 124    2 486.96 IMG_170906_182519_0032   2 0.1684053
## 125   18 305.90 IMG_170906_182519_0032   2 0.1514488
## 126   33 151.64 IMG_170906_182519_0032   2 0.2076212
## 127    2 486.96 IMG_170906_182519_0032   3 0.2052642
## 128   18 305.90 IMG_170906_182519_0032   3 0.2016440
## 129   33 151.64 IMG_170906_182519_0032   3 0.2221849
## 130    2 486.96 IMG_170906_182519_0032   4 0.2000579
## 131   18 305.90 IMG_170906_182519_0032   4 0.1999472
## 132   33 151.64 IMG_170906_182519_0032   4 0.2175644
## 133   18 305.90 IMG_170906_182525_0033   1 0.3112594
## 134   33 151.64 IMG_170906_182525_0033   1 0.5144554
## 135   18 305.90 IMG_170906_182525_0033   2 0.1594713
## 136   33 151.64 IMG_170906_182525_0033   2 0.2494634
## 137   18 305.90 IMG_170906_182525_0033   3 0.2018453
## 138   33 151.64 IMG_170906_182525_0033   3 0.2550449
## 139   18 305.90 IMG_170906_182525_0033   4 0.2235672
## 140   33 151.64 IMG_170906_182525_0033   4 0.2828100
## 141    3 184.99 IMG_170906_182531_0034   1 0.3149959
## 142   19 292.94 IMG_170906_182531_0034   1 0.3125520
## 143   33 151.64 IMG_170906_182531_0034   1 0.5087944
## 144    3 184.99 IMG_170906_182531_0034   2 0.1513753
## 145   19 292.94 IMG_170906_182531_0034   2 0.1535808
## 146   33 151.64 IMG_170906_182531_0034   2 0.2157703
## 147    3 184.99 IMG_170906_182531_0034   3 0.1921439
## 148   19 292.94 IMG_170906_182531_0034   3 0.2061476
## 149   33 151.64 IMG_170906_182531_0034   3 0.2267226
## 150    3 184.99 IMG_170906_182531_0034   4 0.1936413
## 151   19 292.94 IMG_170906_182531_0034   4 0.2079446
## 152   33 151.64 IMG_170906_182531_0034   4 0.2290599

Scatterplots of band values (din) with measured chlorophyll.

toplo <- exts %>% 
  group_by(site, chl, bnd) %>% 
  summarise(din = mean(din))
ggplot(toplo, aes(x = din, y = chl)) + 
  geom_point() + 
  facet_wrap(~bnd, ncol = 2) + 
  stat_smooth(method = 'lm') + 
  theme_bw()